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C*~) The masses of the supermassive black-holes found in galaxy bulges 
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are correlated with a multitude of galaxy properties^ leading 
to suggestions that galaxies and black-holes may evolve together 3 . 
The number of reliably measured black-hole masses is small, and 
the number of methods for measuring them is limited 4 , holding 
back attempts to understand this co-evolution. Directly measur- 
ing black-hole masses is currently possible with stellar kinemat- 
ics (in early-type galaxies), ionised-gas kinematics (in some spiral 
and early-type galaxies 5 7 ) and in rare objects which have central 
maser emission 8 . Here we report that by modelling the effect of 
a black-hole on the kinematics of molecular gas it is possible to 
fit interferometric observations of CO emission and thereby accu- 
rately estimate black hole masses. We study the dynamics of the 
gas in the early-type galaxy NGC 4526, and obtain a best fit which 
requires the presence of a central dark-object of 4.5+g f, x 10 s M 
(3(j confidence limit). With next generation mm-interferometers 
(e.g. ALMA) these observations could be reproduced in galaxies 
out to 75 megaparsecs in less the 5 hours of observing time. The 
use of molecular gas as a kinematic tracer should thus allow one 
to estimate black-hole masses in hundreds of galaxies in the local 
universe, many more than accessible with current techniques. 

For this study, we targeted the fast-rotating early-type galaxy 
(ETG) NGC 4526. This object has a stellar velocity dispersion within 
one effective radius^ o e = 222 km s _1 . This galaxy has not had its 
black-hole mass measured by any other technique to date but, given 
that most galaxies with bulges seem to have a black hole 2 ^, from the 
er - Mbh relation^ we can estimate M BH ~ 2 x 10 s A/ Q . The SMBH 
sphere of influence (SOI) would thus be rsoi = GMbh/c 2 ~ 20 pc 
or C//25 at the 16.4 Mpc distance of NGC 4526 (derived from surface- 
brightness fluctuations"). NGC 4526 also has a large molecular gas 
reservoir (total molecular hydrogen mass Mh, = 3.69 x 10 8 Mgf^j), 
previously mapped at low spatial resolution^ (~4"). The molecular 
gas co-rotates with the stars 14 and Hubble Space Telescope (HST) im- 
ages show that it is coincident with regular dust lanes extending up to 
the galaxy centre, strongly suggesting that molecular gas with regular 
kinematics exists around the SMBH 7 . 

D 

NGC 4526 was observed in the CO(2-l) line (230 GHz) using the 
Combined Array for Research in Millimetre Astronomy (CARMA) in 
A, B a nd C configurations 15 . The data were reduced in the standard 
manne r 6 ' 17 !, ant ] more details are presented in Figure 1. The spatial 
resolution achieved along the kinematic major axis of the galaxy is 
0'.'25 (20 pc), equal to the predicted SMBH SOI. The gas in this source 
seems to be distributed in a central component, and an inner ring with 
spiral spurs which lead outwards to another ring at larger radii. A de- 
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tailed analysis of the molecular gas morphology will be conducted in a 
future work. We note, however, that the gas structures lie in the plane 
of the galaxy and appear to be regularly rotating; they are therefore 
likely to be dynamically cold, and a good tracer of the potential. 

We created a grid of simulations of NGC 4526 spanning a range of 
SMBH masses (no SMBH, Mbh= 3x10 6 , and then 5xl0 7 - 1.45 x 
10 9 Mq in linear steps) and J-band mass-to-light ratios [M/Li = 
0.55 - 6.15 Mq/Lq in linear steps). For full details of these simula- 
tions see Section 1.1 in the online supplementary information. We fix 
the inclination of the gas dislP2l(i — 79°), and use an axisymmetric 
mass model of NGC 4526 9 (carefully fitted to avoid contamination due 
to dust; see Section 1.1.2 in the supplementary material) to derive the 
circular-velocity curve expected from the luminous matter alone. The 
presence of a SMBH in NGC 4526 manifests itself as an inner Kep- 
lerian rise of the rotation curve (above that expected from luminous 
matter only). On larger angular scales such fast-rising rotation curves 
have been observed, and used to infer the masses of central star clus- 
ters and bulges 20 . We fit these models to our observed data in order to 
determine if such an excess due to a central dark mass is detectable in 
NGC 4526. 

Figure [T] shows three different simulated position-velocity dia- 
grams (PVDs) overlaid on the observed PVD of NGC 4526, with 
±1'.'15 insets. The first panel shows the best model with no SMBH, 
with a clear excess of high-velocity molecular gas at the centre. The 
middle panel shows our overall best-fit model, clearly reproducing bet- 
ter the observed PVD at all radii, with Mbh = 4.5 x 10 s Mq and 
M/Lj — 2.65 Mq/Lq. The final panel shows a model with a much 
larger SMBH, clearly incompatible with our data. The lower panels of 
Figure[T]show the trace extracted from the observed PVD (as described 
in the supplementary material), and associated errors. The grey line 
shows the trace extracted from the models in the same way (with points 
denoting the values at the same radius as the observed points). The 
residuals between the data and the model at each position (AV=data 
minus model) are shown in the bottom panel. Clearly the best-fit model 
produces a significantly better fit to the data. 

Figure [2] shows the \ 2 contours of our fits (conducted as de- 
scribed in the supplementary material) as a function of A/bh and 
Ml Li. A clear global minimum is present at the best-fit values, with 
a minimum reduced \ 2 of 1.27. We define uncertainties for each 
of our fitted parameters from likelihood functions of each parame- 
ter (marginalised over the other), by finding the region containing re- 
spectively 68% and 99% of the probability. Our final best-fit values 
are M BH = 4.5tJ'j x 10 8 Mq with M/Lj = 2.65±0.21 Mq/Lq 
(68% confidence level), and 4.5t£? x 10 8 Mq with M/Li = 
2 -65±o:52 Mq/Lq (99% confidence level). These values are fully 
consistent within the 3<r-uncertainties with the reported <j c -Mbh 
relation 10 . Our formal uncertainties (0.14 dex and 0.5 dex at lcr and 
3a, respectively) are also similar (and in many cases smaller) to the av- 
erage uncertainties reported by other authors. For instance the mean lcr 
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Figure 1 | NGC 4526 kinematic models and data. Top: Model PVDs (black contours), overlaid on the observed CO(2-l) PVD (filled orange contours). This PVD 
was created from our CO(2-l) observations of NGC 4526 from CARMA. The synthesized beam size achieved in these observations is 0'.'27 X 0'.'17, and the velocity 
channel width is 10 km s _1 . The final fully reduced and calibrated data cube has an RMS noise of 2 . 88 mJy beam ~ 1 . The PVD was created by rotating the data cube 
to align the kinematic major axis of the molecular gas with the x-axis, and then summing over one beam width around the axis in the y direction. The spatial resolution 
achieved in the PVD is 0'.'25 (20 pc), equal to the predicted SMBH SOI. Our results do not depend on the method used to extract the PVD. From left-to-right, the best 
model with no SMBH, the overall best-fit model, and a model with an overweight SMBH. The model Mbh and M / Lj are indicated in the top-left corner of each 
panel, and an inset of the central ±l"l5 is shown in the bottom-right corner. Middle: Black points show the trace extracted from the observed PVD, and associated 
standard errors. The grey-line shows the trace extracted from the models (and the grey crosses denote the value of the trace at the same radius as the observed points). 
Bottom: Residuals between the model and data at each position (AV=data minus model; km s -1 ). The error bars shown in the middle and bottom plots correspond to 
the formal uncertainties in fitting the trace (see Section 1.2 in the supplementary information), added in quadrature with two factors of 5 km s _1 (to account for the 
finite velocity channel width in both the data and model). 



error reported when using dynamical stellar and ionised gas techniques 
is «0.6 dex-^ (however the rather different systematic errors involved 
make direct comparison difficult). 

The x 2 contours show the usual degeneracy between SMBH mass 
and mass-to-light ratio found in other studies. To fit the data without a 
SMBH, a negative M/L gradient would be necessary, with older stellar 
populations dominating in the inner parts (high M/L) and young stars 
at large radii (low M/L). While the real M/L is unlikely to be exactly 
constant, this is opposite to the trend reported by reddening free stellar 
population studies" 1 and that expected from the presence of molecular 
gas and star-formation in the central regionP^. 

The ability to determine a black-hole mass accurately using 
molecular-gas can be affected by many of the same issues that af- 
fect measurements of ionised-gas. Turbulent motions are in general 
small in molecular gas, but could conceivably increase around a black- 
hole. Similarly, if the inner gas were to be misaligned from the stel- 
lar body, our mass estimate would be systematically affected in a way 
that is degenerate with a change in M/L. In this galaxy, however, 
we find no evidence that the velocity dispersion increases in the inner 
regions (as described in Section 1.1.1 in the supplementary material), 
and constrain the inner gas to be aligned with the stellar body within 
<3°. Such a misalignment could change velocities by an insignificant 
amount ( < 3km s _1 ; see Section 1.1.2 in the supplementary material). 
The presence of dust could also cause mass models to underestimate 
the contribution of luminous matter to the potential. We again believe 
this should not introduce significant errors in this object, because of 
our careful treatment of dust in the mass model (see Section 1.1.2 in 
the Supplementary Information). Future studies using this technique 
should choose their targets to minimise the impact of such effects, but 
can also include warps and turbulent motions in their gas disk models 23 , 



and use near-infrared photometry as done in previous studies of dusty 
objects 2 ^! 

The use of molecular gas as a kinematic tracer holds great promise 
to increase the total number of SMBH mass measurements in galax- 
ies of all types where CO is detected. Angular resolutions of up to 
0'.'15 can be achieved with current mm-interferometers, allowing for 
instance to resolve the SOI of a 2 x 10 8 M Q SMBH in a galaxy with 
<7 e = 200 km s~ only out to around 30 Mpc. The next generation of 
millimetre interferometers will have over an order of magnitude greater 
angular resolution (for example, « 6 mas at 690 GHz with ALMA), as 
well as greatly increased sensitivity. In less the 5 hours of integra- 
tion time with ALMA one could achieve the same sensitivity and lin- 
ear resolution as the observations presented here, in a galaxy 75 Mpc 
away. Galaxies of lower mass (with smaller SMBHs) will also be- 
come accessible. For instance, the SOI of a Milky Way-like SMBH 24 
(M BH = 4 x 10 6 M , cr c = 105 km s" 1 ) will be resolved up to 
« 50 Mpc. 

Measurements using a single technique will thus be possible over 
the entire range of the ct c -A/bh relation, leading to much needed re- 
duced uncertainties in the slope and normalisation of black hole-galaxy 
relations. This will allow the study of a large number of spiral galax- 
ies, that cannot currently be probed easily using stellar dynamical tech- 
niques. Furthermore, it will also provide access to a larger number of 
systems than ionised gas techniques (currently limited by the ~0'.'05 
resolution of HST). Even considering that only some of the accessi- 
ble objects will have suitably relaxed and centrally-peaked distribu- 
tions of molecular gas, measuring SMBH masses will be possible in 
many hundreds of spiral (and early-type) galaxies, many times more 
than possible using conventional stellar and ionised-gas tracers today. 
This technique could also be extended using other spectral lines, such 
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Figure 2 | NGC 4526 SMBH mass uncertainties. Ax 2 = x 2 - Xmlnimum 

contours of our fits to the CO(2-l) PVD, as a function of the two free parameters 
Mbh and M/Lj. The model grid is shown with black dots and the overall best- 
fit model by a white star. The solid shaded contours with black lines correspond 
to the l-3cr levels with one degree-of-freedom (Ax 2 = 1, 4, 9) with good models 
in the darkest areas. The 4 & 5<r levels (A x 2 = 16, 25) are denoted with dashed 
grey lines. 



as the very bright atomic gas cooling lines at higher frequencies (in 
order to increase the spatial resolution achieved even further), and opti- 
cally thin transitions such as 13 CO (which would allow determinations 
of black-hole masses even in systems that are exactly edge on). 
Received 11th September; accepted 27th November 2012 
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Supplementary Information 

1.1 Model Observations 

In this paper we compare the observed CO position- 
velocity diagram with one created by modelling the 
gas disk. When creating these models we included the 
observational effects of beam-smearing, disk-thickness, 
velocity dispersion, binning, etc, by utilising the KINe- 
matic Molecular Simulation (KinMS) mm-wave obser- 
vation simulation tool 16 ' 25 . The simulations used a 
beam, pixel size and velocity resolution identical to our 
observations, and an idealised surface brightness profile 
which fits the data. We use an axisymmetric mass model 
of NGC 4526 to derive the circular-velocity curve ex- 
pected from the luminous matter alone. These elements 
of the models are described in detail below. 

1.1.1 Gas distribution 

The distribution of the molecular gas in this source is 
a required input for the KinMS tool. As we primarily 
aim to assess the kinematic agreement between the data 
and models with different SMBH masses (rather than 
the details of the gas distribution itself), a simple sur- 
face brightness profile fitting the observed distribution 
is sufficient. We therefore adopted a profile composed 
of a central Gaussian with full-width-at-half-maximum 
FWHM=0"15« 12 pc, centred at R=0.6"« 48 pc (nu- 
clear ring) and a Gaussian of FWHM=2"« 160 pc cen- 
tred at R=5"« 400 pc (main ring). This parameterized 
surface brightness profile fits well the radial distribution 
of the observed clean components within the inner 6", 
and identical results are obtained for all sensible profile 
choices. We assumed a gas velocity dispersion of 8 km 
s~ 1 here (our results do not depend on choosing any rea- 
sonable value for this parameter). Models with er gas < 
10 km s _1 are required to match the observed molecular 
gas distribution at all radii. 

1.1.2 Mass Model 

The mass model we use in this work is based on a 
multi-Gaussian expansion (MGE 26 ) of F8 14W HST and 
MDM 1.3-m telescope /-band images (probing both 
small and large scales), and is parametrised by the I- 
band stellar mass-to-light ratio ftl/Li and inclination i. 
From it the circular velocity of the galaxy can be com- 
puted at any radius. The model is carefully fitted to re- 
move any contamination from dust still visible in the 
/-band. The dust disk in this system is well defined, and 
we have many resolution elements inside the SOI of the 
SMBH. This allows us to carefully mask the affected re- 
gions when conducting the MGE fit 27 . The unaffected 
sectors were then used to derive the MGE model. Such 
a process has been shown to work well at recovering the 
intrinsic light distribution in such systems 6 . 



Axisymmetric two-integral Jeans and three-integral 
Schwarzschild dynamical models of the large-scale 
stellar kinematics (within 45") suggest M/Lj = 
3.35 Mq/Lq and i = 79°, but stellar population 
models of the inner parts (^5") suggest M/Lj = 
2.6 Mq/Lq (derived assuming a Kroupa IMF 6 , which 
is likely appropriate for this galaxy 28 ), presumably be- 
cause of young stellar populations associated with the 
molecular gas. Any residual dust obscuration in the very 
centre of the galaxy could cause us to underestimate the 
M/L in the dynamical model, however the stellar popu- 
lation Mj Lj estimate (derived from an analysis of line- 
strength indices close together in wavelength 18 ) is dust- 
independent. 

The inclination estimated from the dynamical models 
is consistent with that derived for the dust disk and the 
spread of measured values is very small 17 (±3°). 

Ellipse fitting to the inner dust lanes visible in HST 
imaging, and to the inner molecular component show 
that the inner gas is also aligned with the stellar body 
within <3° As NGC 4526 is close to edge on, such an 
uncertainty in inclination produces an insignificant er- 
ror in the projected velocities (w 3 km s" 1 ), which is 
insignificant when compared to the uncertainties intro- 
duced by other parameters (see Section 1.2 of the sup- 
plementary material). We therefore fix the inclination of 
the gas disk at all radii to i = 79°. 



1.2 Trace extraction 

For a quantitative comparison focusing on the kinemat- 
ics, we extracted the 'trace' of the outer envelope of the 
observed and simulated PVDs, every O'.'l in radius (in 
order to Nyquist sample our beam), by fitting a Gaussian 
(along the velocity axis) to each spatial bin and taking 
the central velocity of the Gaussian plus one standard 
deviation as the trace. The uncertainties were estimated 
from the formal uncertainties in the velocity and stan- 
dard deviation of the fitted Gaussian, plus two factors 
of 5 km s _1 (to account for the finite velocity channel 
width in both the data and model). These channel width 
factors usually dominate. 



1.3 Fitting procedure 

The best fit was found using a standard \ 2 test, compar- 
ing the one-dimensional trace of the data and simula- 
tions over the central ±4'/5, the area most affected by a 
putative SMBH (See Figure 2). We did not fit inside the 
innermost beam, where the effects of beam smearing are 
strongest, leaving us with 68 observed data-points. The 
best-fit profile has a x 2 =86.36, equivalent to a reduced 
X 2 =1.27. 
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